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Coalescing binaries of neutron stars and black holes are one of the most important sources of gravitational 
waves for the upcoming network of ground based detectors. Detection and extraction of astrophysical infor¬ 
mation from gravitational-wave signals requires accurate waveform models. The Effective-One-Body and other 
phenomenological models interpolate between analytic results and numerical relativity simulations, that typi¬ 
cally span 0(10) orbits before coalescence. In this paper we study the faithfulness of these models for neutron 
star - black hole binaries. We investigate their accuracy using new NR simulations that span 36 — 88 orbits, 
with mass-ratios q and black hole spins xbh of ( q,XBH ) = (7, ±0.4), (7, ±0.6), and (5, —0.9). We find 
that: (i) the recently published SEOBNRvl and SEOBNRv2 models of the Effective-One-Body family dis¬ 
agree with each other (mismatches of a few percent) for black hole spins xbh > 0.5 or xbh < —0.3, with 
waveform mismatch accumulating during early inspiral; (ii) comparison with numerical waveforms indicate 
that this disagreement is due to phasing errors of SEOBNRvl, with SEOBNRv2 in good agreement with all of 
our simulations; (iii) Phenomenological waveforms agree with SEOBNRv2 only for comparable-mass low-spin 
binaries, with overlaps below 0.7 elsewhere in the neutron star - black hole binary parameter space; (iv) com¬ 
parison with numerical waveforms shows that most of this model’s dephasing accumulates near the frequency 
interval where it switches to a phenomenological phasing prescription; and finally (v) both SEOBNR and post- 
Newtonian models are effectual for neutron star - black hole systems, but post-Newtonian waveforms will give 
a significant bias in parameter recovery. Our results suggest that future gravitational-wave detection searches 
and parameter estimation efforts would benefit from using SEOBNRv2 waveform templates when focused on 
neutron star - black hole systems with q < 7 and xbh ~ [—0.9, ±0.6]. For larger black hole spins and/or 
binary mass-ratios, we recommend the models be further investigated as NR simulations in that region of the 
parameter space become available. 


I. INTRODUCTION 

The Advanced Laser Interferometer Gravitational-wave 
Observatory (aLIGO) [1, 2] is currently being commissioned 
and will begin observation in 2015, reaching its design sensi¬ 
tivity by 2018 — 19. The Virgo gravitational-wave observa¬ 
tory [3] will begin operation in 2016. With improved sensi¬ 
tivity, these detectors will access a thousand times as much 
volume as their first generation counterparts. In addition, the 
KAGRA detector is currently under construction in Japan [4], 
and a plan to build an advanced LIGO detector in India is un¬ 
der consideration. Compact binaries are the most promising 
sources of gravitational waves (GWs) for aLIGO. Binary sys¬ 
tems containing stellar-mass black holes (BH) and/or neutron 
stars (NS) inspiral and merge because of their GW emission. 
The GW waves emitted with frequencies above ~ 10 Hz will 
be in the sensitive band of aLIGO and Virgo. 

In this paper, we focus on neutron star - black hole (NSBH) 
binaries. Based on our current understanding of the astro- 
physical NS and BH population, stellar binary evolution, and 
on population synthesis studies, we expect aLIGO to observe 
0.2 — 300 NSBH binary mergers a year [5]. GW observations 
of NSBH binaries have significant scientific potential, beyond 
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the initial discovery of a new class of astrophysical systems. 
GWs emitted by coalescing NSBH binaries carry signatures 
of strong-field gravitational dynamics. Unlike binary neutron 
stars, GWs from NSBH binaries will contain the signatures 
of the interaction of BH spins [6-11] with the orbital mo¬ 
tion. Significant efforts are underway to access this infor¬ 
mation using the aLIGO and Virgo detector network to test 
general relativity in the strong gravity regime [12, 13]. The 
observation and characterization of a population of NSBH 
sources will also shed light on stellar evolution and compact¬ 
binary formation mechanisms: e.g., a gap in the mass distri¬ 
bution of NSs and BHs could shed light on the mechanism 
of supernova explosions [14-16]. An unambiguous detec¬ 
tion of GWs from a NSBH system accompanied by electro¬ 
magnetic observations could provide information about the 
internal structure of NSs [17] and could provide strong evi¬ 
dence linking compact binary mergers and short Gamma-ray 
bursts (SGRBs) [18-21]. However, unlocking the full scien¬ 
tific potential of GWs emitted by NSBH coalescences requires 
both detecting as many of such signals as possible and accu¬ 
rately characterizing them to understand the properties of their 
source binaries. 

Detection searches are based on the matched-filtering tech¬ 
nique [22], using modeled waveforms as filter templates. 
Searches for compact binaries with initial LIGO and Virgo de¬ 
tectors used non-spinning template waveforms [23-27] (with 
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the exception of [28], ChrisVDBrock et al). While the cata¬ 
loged astrophysical population of NSs have small spins (mass- 
normalized |x| < 0.05), the spins of stellar-mass BHs are un¬ 
certain, with estimates ranging from low to nearly extremal 
values (i.e., nearly as fast as possible—see, e.g. [11, 29, 30] 
for examples of nearly extremal estimates of BH spins, and 
see Refs. [31, 32] for recent reviews of astrophysical BH spin 
measurements). Recent work has shown that including non- 
precessing (that is, aligned) component spins in templates 
used in matched-filtering gravitational-wave searches will sig¬ 
nificantly improve the searches’ sensitivity [33]. Therefore, 
aLIGO-Virgo searches targeting NSBH binaries plan to use 
aligned-spin waveform templates [34]. Because they are cen¬ 
tral to matched-filtering searches, it is crucial to have GW 
models that accurately capture the NSBH coalescence pro¬ 
cess. Modeling inaccuracy would reduce the signal-to-noise 
ratio (SNR) recovered by detection searches, and degrade the 
range of aLIGO-Virgo observatories. It would also lead to 
systematic, but not necessarily controlled, biases in the recov¬ 
ered masses and spins of the source. 

Studies of the accuracy of contemporary waveform mod¬ 
els in the past have focused on post-Newtonian (PN) [35] 
and recent Effective-One-Body (in particular, the “SEOB- 
NRvl” [36]) models. It has been shown that PN approxi- 
mants disagree significantly with each other and with SEOB- 
NRvl for aligned-spin NSBH binaries [37], despite the in¬ 
clusion of the highest-known order spin contributions to the 
binary phasing [38, 39]. While the accuracy of the SEOBNR 
models is enhanced through calibration against high-accuracy 
Numerical Relativity (NR) merger simulations, most of these 
simulations correspond to comparable mass ratios. Therefore, 
the extension of SEOBNR into NSBH parameter space is not 
guaranteed to be reliable. 

In this paper, we systematically investigate waveform ap- 
proximants in the context of NSBH binaries. Unlike past 
studies, we investigate not just the precision (mutual agree¬ 
ment of approximants) but also their accuracy , by comparing 
with long NR simulations with q = mBu/m^s = {5, 7} and 

aligned BH spin xbh = 5bh/-M| h = {±0.4, ± 0 . 6 ,-0.9}. 

(Note that, except where we specify otherwise, we adopt ge- 
ometrized units with G = c = 1 in this paper). These sim¬ 
ulations are described further in Sec. II. In addition to PN 
and SEOBNRvl, we compare with the more recent SEOB¬ 
NR v2 and the phenomenological PhenomC [40] models. We 
use the zero-detuning high power noise curve for Advanced 
LIGO [41] with a 15 Hz lower frequency cutoff in our calcula¬ 
tions. We allow the BH spin to vary over [—1,1], and its mass 
to vary over [3M Q , 14M 0 ]. The NS mass is fixed at ttins = 
1.4M e with xns = 0, as is consistent with the observed astro- 
physical NS population [42^14]. Note that while investigating 
waveform modeling errors, we ignore NS matter effects and 
treat the NS as a low-mass BH. Although matter effects are 
expected to be measurable by aLIGO (e.g. [45, 46]), they af¬ 
fect the inspiral phasing starting at 5PN order. As there are 
lower-order spin-dependent vacuum terms in PN phasing that 
remain unknown, the effect of ignoring matter-dependent sec¬ 
ular terms will be sub-dominant to other sources of error in 
waveform models. We also ignore the effect of NS disruption 


before merger, which is likely when the mass-ratio rriBu/m^s 
is small and/or the BH has relatively high aligned spin [47]. 
However, this disruption occurs at fairly high frequencies, i.e. 
at /gw > 1-2 kHz [47], and its effects are expected to be 
small due to the significantly reduced sensitivity of aLIGO at 
such frequencies [48]. We leave the study of this effect to 
future work. 

First, we study GW model precision by comparing the 
PN time-domain TaylorT4, PN frequency-domain TaylorF2, 
SEOBNRvl and PhenomC models with the most recent 
SEOBNRv2 model. As SEOBNRv2 has been calibrated to 38 
NR simulations, we take it as the fiducial model representing 
the true waveform. We find that both PN models have overlaps 
with SEOBNRv2 below 0.9 for mass-ratio q > 3 and/or BH 
spin | xbh | > 0.5. We also find that PhenomC and SEOB¬ 
NR v2 have overlaps below 0.9 for q < 5 and/or BH spin 
|Xbh| > 0-3, falling as low as 0.6. Finally, we also find the 
overlaps between SEOBNRvl and SEOBNRv2 fall below 0.9 
for NSBH systems with anti-aligned BH spins xbh < —0.5. 

We further investigate the accumulation of mismatch be¬ 
tween different models, as a function of GW frequency. For 
PN approximants, we find that most of the mismatch is ac¬ 
crued during the late-inspiral phase when the PN velocity pa¬ 
rameter v/c > 0.2. This is expected, because PN results are 
perturbative expansions in v/c that break down when v be¬ 
comes comparable to c near the time of coalescence. Between 
SEOBNRvl and SEOBNRv2, we find that mismatch is ac¬ 
crued during the early-to-late inspiral transition period when 
v/c < 0.26. We find a similar trend between PhenomC and 
SEOBNRv2. This demonstrates discrepancies between NR- 
calibrated models in the early inspiral phase, despite good 
agreement close to merger, where all of the models have been 
calibrated to NR. 

Second, we study the accuracy of NSBH waveform models 
by computing their overlaps (or faithfulness) against our long 
NR simulations. Our simulations extend down to v/c ~ 0.2, 
and are long enough to probe the frequency range in which 
SEOBNRvl/PhenomC phase evolutions differ from SEOB- 
NRv2. We find that SEOBNRvl has 1 — 3% mismatches 
against the aligned-spin simulations, which rise up to ^ 5% 
against the anti-aligned-spin ones. While most of this mis¬ 
match is accumulated during the last few pre-merger orbits 
for aligned -spin cases, for anti-aligned cases it accumulates 
over the 30 — 50 inspiral orbits that our simulations span. On 
the other hand, we found that SEOBNRv2 has < 1% mis¬ 
matches with NR, for both aligned and anti-aligned simula¬ 
tions. Therefore, we conclude that the differences between 
the two SEOBNR models are because of the phasing errors 
in SEOBNRvl. For PhenomC and both PN approximants, 
we find > 10% mismatches against NR, for both aligned 
and anti-aligned spin simulations. We therefore conclude 
that SEOBNRv2 provides the most accurate description of 
aligned-spin NS-BH coalescence waveforms, with the caveat 
that the model should be analyzed for more extreme compo¬ 
nent spins. 

Third, we investigate the suitability of different models for 
detection searches. To address this question, we compute 
the effectualness of different models by allowing the addi- 
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tional degree of freedom of maximizing overlaps between an¬ 
alytic and numerical waveforms over intrinsic binary parame¬ 
ters. We find that both SEOBNRvl and SEOBNRv2 recover 
> 99.8% of the optimal SNR. PhenomC shows low SNR re¬ 
covery, which drops below ~ 90% for anti-aligned BH spins. 
We therefore recommend against using this model in NSBH 
detection searches. Both PN models recover about 98% of the 
SNR for aligned-spin systems, and are therefore effectual for 
aLIGO searches. For anti-aligned systems, both TaylorT4 and 
TaylorF2 models recover < 96% of the SNR and would likely 
benefit from the computation of higher order spin-dependent 
corrections to PN dynamics. Therefore we recommend that 
SEOBNRv2 be preferred in aLIGO NSBH detection searches. 

Finally, we probe the question of systematic biases 
in parameter recovery corresponding to using each ap- 
proximant to model aLIGO parameter estimation tem¬ 
plates. We find that the accuracy of the chirp mass 
(A4 C = (raira 2 ) 3 / 5 (rai + m 2 ) _ly/5 ) recovery increases with 
the number of orbits that are integrated over. All approxi- 
mants recovered A4 C within a few percent. The spin-mass- 
ratio degeneracy makes the accurate determination of mass- 
ratio and component spins more challenging. We find sys¬ 
tematic biases in mass-ratio to be between a few to tens of 
percents, increasing with BH spin, with similar biases in the 
recovered values of BH spins. Of all the models considered, 
we find that SEOBNRv2 surpasses others in faithfulness. Us¬ 
ing the accuracy measures proposed in [49], we also found 
SEOBNRv2 to be indistinguishable from true waveforms up 
to SNRs « 8 — 14(16 — 18) for aligned (anti-aligned) BH 
spins. We therefore recommend that SEOBNRv2 be used in 
aLIGO parameter estimation efforts for aligned-spin NSBH 
detection candidates, but we also recommend that SEOBNR 
be tested for higher component spins. 

Our results are limited by the fact that our NR waveforms 
only extend down to v/c ~ 0.21 — 0.24 (i.e. 60 — 80 Hz 
for NSBH masses), while aLIGO is sensitive down to 15 Hz. 
A sizable fraction (35 — 45%, depending on BH spin) of the 
signal power will be accumulated at frequencies below this 
range. We plan to extend these calculations to lower frequen¬ 
cies in future work. 

The remainder of the paper is organized as follows: In 
Sec. II, we describe the NR waveforms presented in this pa¬ 
per and discuss their convergence. In Sec. Ill, we describe 
the waveform models studied here. In Sec. IV, we describe 
the measures used to quantify waveform discrepancies. In 
Sec. V, we discuss the faithfulness of different waveform ap- 
proximants for different NSBH masses and spins, and also as 
a function of the emitted GW frequency. In Sec. VI, we inves¬ 
tigate the late-inspiral accuracy of all approximants using our 
high accuracy numerical simulations. In Sec. VII, we study 
the viability of using different approximants as detection tem¬ 
plates, as well as their intrinsic parameter biases for aLIGO 
parameter estimation studies. In Sec. VIII, we summarize and 
discuss our results. 


II. NUMERICAL RELATIVITY SIMULATIONS 

We construct our NR waveforms using the Spectral Ein¬ 
stein Code (SpEC) [50]. Their parameters are summarized 
in Table I. Of particular note is the length of these simula¬ 
tions; the shortest waveform presented here has over 36 orbits 
of inspiral before the merger and the longest has nearly 90 or¬ 
bits. Fig. 1 shows the real part of the waveform dimenionless 
strain, rh 22 /M for each of the simulations on Table I. With 
the exception of the 176 orbit simulation presented in [51] 
and a 48.5 orbit simulation presented in [52] these waveforms 
are among the longest done to date. The longest waveform 
currently in the SXS catalog is only 35.5 orbits [53]. 

To test the accuracy of the simulations we ran each simula¬ 
tion using different numerical resolutions; we label each reso¬ 
lution by an integer AT, where larger AT indicates finer resolu¬ 
tion. We compute the phase of the £ = 2, m = 2 mode of 4/ 4 
(the second time derivative of the complex strain /i 22 ^)) for 
different resolutions Fig. 2 shows these phase differences for 
each pair of resolutions for five of the seven numerical simu¬ 
lations presented here. (We ran the other two simulations at 
fewer than 3 values of AT, so such a plot would not be useful 
for those cases.) 

If for all subdomains, i) the number of grid points increased 
uniformly with increasing AT, and ii) at any given time, the 
locations of the boundaries of all subdomains were indepen¬ 
dent of AT, then Fig. 2 would represent a classic convergence 
test. In that case, the phase differences should decrease with 
AT in a predictable way, according to the convergence order 
of the numerical scheme. The simulations here, however, use 
a spectral adaptive mesh refinement (AMR) scheme [54], and 
the label AT determines the error tolerance used by AMR when 
it decides whether to change the number of points in a given 
subdomain and when it decides whether to split a single sub- 
domain into many smaller ones or to join several subdomains 
into a larger one. Because AMR makes these decisions in¬ 
dependently for different values of AT, at any given time it is 
possible that a given subdomain has the same number of grid 
points for two values of AT, and it is possible that subdomain 
boundaries for different values of AT do not agree. There¬ 
fore, we do not necessarily expect strict convergence in Fig. 2. 
These issues will be discussed in detail in a separate paper that 
focuses on convergence of BBH runs using SpEC. 

Nevertheless, for the x = ±0.6 simulations, the differ¬ 
ences converge well with AT : differences become successively 
smaller with increasing resolution. Furthermore, the differ¬ 
ence between AT = 3 and TV = 2 is approximately equal to 
the difference between N = 4 and AT = 2, indicating that 
these differences essentially measure the error in N = 2. For 
the x = ±0.4 and x = —0.9 simulations, the difference be¬ 
tween AT = 3 and AT = 2 is smaller than the difference be¬ 
tween N = 3 and AT = 1, but the spacing between differences 
is not uniform, and there are some anomalously small phase 
differences, such as between AT = 2 and AT = 1 for x = 0.4. 

For x = ±0.6, the difference between the two finest resolu¬ 
tions AT = 3 and AT = 4 is a good measure of the numerical er¬ 
ror in the AT = 3 simulation. The error in AT = 4 could be sim¬ 
ilarly measured via the difference between AT = 4 and AT = 5, 
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ID 

q 

Xi 

No. of orbits 

Initial f gw (Hz) 

Initial Mu or bitai 

a 

D 0 

e 

SXS:BBH:0202 

7 

(0,0,0.6) 

62.1 

75.5 

0.01309 

4.3970 x 10“ b 

17.0000 

9 x 10“ b 

SXS:BBH:0203 

7 

(0,0,0.4) 

58.5 

76.0 

0.01317 

-9.8403 x 10“ 6 

17.0005 

< 1.6 x 10“ 4 

SXS:BBH:0204 

7 

(0,0,0.4) 

88.4 

60.3 

0.01045 

-4.6373 x 10“ 6 

20.0000 

< 1.7 x 10“ 4 

SXS:BBH:0205 

7 

jD 

O 

1 

O 

44.9 

76.0 

0.01318 

-1.4760 x 10“ 5 

17.1036 

7.0 x 10 -5 

SXS:BBH:0206 

7 

(0,0, -0.4) 

73.2 

59.8 

0.01036 

-7.8300 x 10“ 6 

20.2167 

< 1.6 x 10“ 4 

SXS:BBH:0207 

7 

(0,0, -0.6) 

36.1 

80.8 

0.01399 

7.1708 x 10“ 6 

16.4000 

1.693 x 10“ 4 

SXS:BBH:0208 

5 

(0,0, -0.9) 

49.9 

80.0 

0.0104 

-4.5088 x 10 -5 

20.0778 

5.074 x 10“ 4 


TABLE I. Numerical-relativity simulations used in this study (each performed using SpEC [50]). For each simulation (labeled by ID), the 
table shows the mass ratio q = mi/m 2 > 1 , the spin xi of the heavier compact object (the lighter object is non-spinning), the number of 
orbits, the initial gravitational-wave frequency f gw when the total mass is scaled so that the system mimics a NSBH binary with a NS mass 
of 1.4M®, the initial dimensionless orbital velocity Muj or bitai, radial velocity d, separation Do, and eccentricity e. These values listed for 
f gw , Muj orbital, cl and Do are for the initial data, before junk radiation. The xi — ±0.4 simulations use CFMS initial data while the rest use 
SKS initial data. 



t/M 


FIG. 1. The real part of rh^/M of the £ — m = 2 mode of the numerical waveforms used in this paper, where M is the total mass 
and r is the radial distance from the source to an observer. The waveform labeled by N correspond to simulation SXS:BBH:7V, where 
N £ {202, 203, 204, 205, 206, 207}. The waves are shown as a function of time t. A constant vertical offset is applied to each waveform for 
clarity, and the waves are offset in time so the peak amplitude occurs at time t = 0. 



FIG. 2. Phase differences between different resolutions for the 
£ = 2,ra = 2 modes of Tg, plotted as a function of time. Only 
the five numerical simulations with more than 2 resolutions are dis¬ 
played. The legends indicate which resolutions are compared, e.g. 
“3-2” compares the phase of N = 3 versus N = 2. 


but since we do not have an N = 5 simulation, we take the 
difference between N = 3 and N = 4 as an extremely con¬ 
servative estimate of the error in N = 4. For x = ±0.4 and 
X = —0.9, where the convergence with N is not so appar¬ 
ent, we likewise take the difference between the two highest 
resolutions as the error estimate of the highest-resolution sim¬ 
ulation; in these cases the error estimate is likely not as con¬ 
servative as for x = ±0.6. 


III. WAVEFORM APPROXIMANTS 

In this paper, we consider three waveform families: 
post-Newtonian, Effective-One-Body, and phenomenological 
models [40, 55]. We briefly summarize them here, pointing 
the reader to the references for more detailed descriptions. 
We consider the (£,m) = (2, ±2) spin-weighted spherical 
harmonic waveform multipoles, since (i) these are the dom¬ 
inant modes for non-precessing systems, with the contribu¬ 
tions from other modes being smaller by a few orders of mag¬ 
nitude, and (ii) none of the contemporary IMR models include 
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the sub-dominant waveform modes. 

Post-Newtonian: The Post-Newtonian (PN) approxima¬ 
tion is a perturbative expansion of compact binary inspiral 
dynamics in the limit of slow motions and weak fields. The 
orbital energy E of a non-precessing binary and the flux F of 
the gravitational energy emitted as GWs are known to 3.5PN 
order [56-64]. Combining E and F using the energy bal¬ 
ance equation dE/dt = — F yields a system of differential 
equations; solving these equations gives the GW phase and 
the orbital frequency evolution. The energy balance equation 
can be re-expanded and solved in different ways to obtain dif¬ 
ferent approximants that agree to 3.5PN order but differ at 
higher orders. The PN formalism and the corresponding equa¬ 
tions of motion break down before merger as the underlying 
approximations (slow motions and weak fields) break down; 
therefore, the PN formalism produces only the inspiral por¬ 
tions of the waveform. In this paper, we will examine two 
particular PN approximants, the time-domain TaylorT4 and 
frequency-domain TaylorF2. In both, we include the recently 
published spin-orbit tail (3PN) and the next-to-next-to-leading 
order spin-orbit (3.5PN) contributions [65, 66 ]. We refer the 
reader to the Appendix of Ref. [37] for a summary of the ex¬ 
pressions that describe both approximants. 

Effective-One-Body: The Effective-One-Body (EOB) ap¬ 
proach solves for the dynamics of a compact binary system by 
mapping them to the dynamics of an effective test particle of 
mass p = 77117712 /(mi -\-rn 2 ) with a spin 5*(mi, m 2 , Si, S 2 ) 
in a space-time described by a suitable deformation of the 
Kerr metric. Specifically, when constructing model wave¬ 
forms using the EOB approach, one chooses the deforma¬ 
tion and the test-particle spin such that the geodesic followed 
by the test particle reproduces the PN-expanded dynamics of 
the compact binary system with component masses mi and 
m 2 and spins Si and S 2 . Then, one matches the coefficients 
of the deformed metric to the PN expansion up to 3PN or¬ 
der; to further improve accuracy, one adds adjustable 4PN 
and 5PN terms that are calibrated by forcing agreement with 
NR waveforms. The conservative dynamics of the test parti¬ 
cle in the deformed-Kerr spacetime are described by the EOB 
Hamiltonian TTeob- The expressions for H^ob for differ¬ 
ent Spinning-EOB (SEOBNR) models differ at high PN or¬ 
ders and can be found in Refs. [36, 67]. Further, a radiation 
reaction term in the equations of motion captures the non¬ 
conservative dynamics, i.e., the motion of the binary through 
inspiral to merger [36, 67]. In contrast to PN waveforms, EOB 
waveforms continue through merger and ringdown, with the 
ringdown waveform constructed as a linear superposition of 
the first eight quasi-normal modes (QNMs) of the Kerr BH 
formed at binary merger [ 68 ]. Matching the ringdown wave¬ 
form to the inspiral-merger portion determines the coefficients 
associated with these QNMs. 

We use two SEOBNR models (both available in the LIGO 
Algorithm Library (LAL) [69]) in this study: SEOBNRvl 
and SEOBNRv2. These models differ in their calibration 
to NR waveforms. SEOBNRvl models the inspiral-merger- 
ringdown (IMR) waveform for binaries with component spins 
—1 < xi ? 2 < 0.6; while SEOBNRv2 can model more ex¬ 
tremal component spins, i.e. — 1 < X < 0.99. SEOB¬ 


NRvl uses 5 non-spinning NR simulations at mass ratios q = 
{1,2, 3,4, 6 } and two equal-mass simulations with xi ,2 = 
±0.4 to choose six adjustable parameters; note that NSBH 
systems are outside the domain of calibration of this model. 
SEOBNRv2, in addition, includes an adjustable parameter in 
the effective particle spin mapping, one in the Hamiltonian, 
and one in the complex phase of ^ 2 , 2 • These parameters have 
been calibrated against 8 non-spinning and 30 aligned-spin 
NR waveforms. We refer the readers to Ref. [36] and Ref. [67] 
for comprehensive descriptions of SEOBNRvl and SEOB- 
NRv2, respectively. 

Phenomenological model: PhenomC is a phenomenolog¬ 
ical model that has a closed form in the frequency domain 
and describes the GW emitted by aligned spin binaries dur¬ 
ing their inspiral-merger-ringdown (IMR) phases [40] . Closed 
form TaylorF2 expressions capture the early adiabatic inspiral 
stage of binary coalescence, with the frequency-domain wave¬ 
form amplitude and phase expressed as expansions in GW 
frequency. PhenomC inspiral phasing includes non-spinning 
terms up to 3.5PN order, spin-orbit and spin-spin coupling 
terms up to 2.5PN order, and horizon absorption terms up to 
2.5PN order [70]. In the late-inspiral/pre-merger regime, the 
PN approximation is insufficient to model the phase evolution 
accurately; instead, a phenomenologically fitted power series 
in frequency (i.e., a polynomial in / 1 / 3 ) captures the phase 
evolution. Calibrating against a set of NR waveforms [40] 
in the frequency range [0.1/rd, /rd], where / RD is the pri¬ 
mary (least damped) GW frequency emitted during the quasi¬ 
normal ringing of the post-merger black hole remnant, deter¬ 
mines the free coefficients in the resulting pre-merger phase 
prescription In the ringdown regime, PhenomC models the 
phase as a linear function in GW frequency, capturing the 
effect of the leading quasi-normal mode. Similarly, Phe¬ 
nomC constructs the amplitude prescription through piece- 
wise modeling of the pre-merger and ringdown regimes, ap¬ 
proximating the amplitude by a power-series in J 1 / 3 pre¬ 
merger and by a Lorentzian post-merger. For a complete de¬ 
scription of PhenomC and its calibration, we refer the reader 
to Ref. [40] . Note that PhenomC belongs to the unique class 
of models that are both closed-form in the frequency domain, 
and include the late-inspiral, merger, and ringdown in the 
waveforms. These features are especially convenient for real 
GW searches, which operate in the frequency domain, fil¬ 
tering observatory data with a large number (lO 5 — 10 6 ) of 
waveform templates. 


IV. QUANTIFYING WAVEFORM ACCURACY 


As a measure of how “close” two waveforms hi and /i 2 
are in the waveform manifold, we use the maximized overlap 
(match) O , defined by 


0(hiM) 


max (/ii|/i 2 (0 c ,£c)) 

0C,£c 




(i) 
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where the overlap (-|-) between two waveforms is 


(hi\h 2 ) = 4 


r 


hiU)KU) 

Snif) 


d/; 


( 2 ) 


(j) c and t c are the phase and time shift differences between hi 
and h 2 \ h(f ) is the Fourier transform of the GW waveform 
h\ S n (f ) is the one-sided power spectral density (PSD) of the 
detector noise, which we assume to be stationary and Gaus¬ 
sian with zero mean; / min is the lower frequency cutoff for 
filtering; and is the Nyquist frequency corresponding to 
the waveform sampling rate. The normalization of O takes 
away the effect of any overall amplitude scaling differences 
between hi and h 2 . The complimentary measure of the mis¬ 
match A4 between the two waveforms is therefore 


M(h u he)*tl~0(h u h2). ( 3 ) 

Matched-filtering based searches use a discrete bank of 
modeled waveforms as filters. The optimal value of the recov¬ 
ered SNR is p 0 pt = where h tr is the actual GW 

signal in the detector output data. With a finite bank of filter 
templates, the recovered SNR is p ~ 0(h tr , h^) p opt < p opt , 
where h^ is the filter template in the bank that has the highest 
maximized overlap with the signal h tr ; i.e. the recovered SNR 
is maximized over both intrinsic (component mass and spin) 
and extrinsic (0 C and t c ) parameters that describe the source 
binary. For a NSBH population uniformly distributed in spa¬ 
tial volume, the detection rate is oc 0(h tr , h^) 3 . To maxi¬ 
mize the detection rate, it is therefore crucial for the model 
waveform manifold, containing to faithfully reproduce the 
manifold of true waveforms that contains h tT . 

In this paper, we take S n (\f\) to be the zero-detuning high 
power noise curve for aLIGO [41, 48]. The peak GW fre¬ 
quency for the lowest binary masses that we consider, i.e. for 
mi + m 2 — 8.4M©, is ~ 3 kHz during ringdown. We sample 
the waveforms at 8192 Hz, preserving the information content 
up to the Nyquist frequency = 4096 Hz. 

Our numerical waveforms begin at relatively low frequen¬ 
cies (between about 60 Hz and 80 Hz), but nevertheless they 
do not completely span the detector’s sensitive frequency 
band. The discontinuity at the start of the NR waveform, be¬ 
cause of Gibbs phenomena [71], corrupts the Fourier trans¬ 
form. We therefore taper the start of the waveforms using a 
cosine tapering window, whose width is chosen to control the 
corruption of the resulting waveform in a way that the mis¬ 
matches because of tapering stay below 0.2%. Additionally, to 
minimize residual spectral leakage, we apply an eighth order 
Butterworth high-pass filter with the cutoff frequency equal to 
the frequency of the waveform at the end of the tapering win¬ 
dow. In Sec. VIB, we measure the effect of waveform condi¬ 
tioning on the NR waveforms by comparing identically con¬ 
ditioned waveforms against unconditioned-but-significantly- 
longer SEOBNRv2 waveforms. 


V. COMPARISON OF WAVEFORM MODELS 

In this section, we show the faithfulness between wave¬ 
forms from different approximants, where we choose the 
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3456789 10 

Mass ratio 

FIG. 3. In this figure, we show the match between the two SEOBNR 
models we consider in this paper, as a function of the BH spin and 
binary mass ratio. The mass of the NS is fixed at 1.4M©, and its spin 
is set to 0, consistent with the observed astrophysical population of 
NSs [42]. The blue (black) solid lines are the contours of 99% (97%) 
match. We find that the matches between the models drop sharply for 
Xbh < —0.15. 



physical parameters to be consistent with NSBH sources. We 
compare the inspiral-merger-ringdown (IMR) models SEOB- 
NRvl and PhenomC, and the PN models TaylorT4 and Tay- 
lorF2, with the SEOBNRv2 model. We also show how the 
disagreement between approximants builds up over the coarse 
of a binary’s inspiral, by computing their faithfulness over dif¬ 
ferent GW frequency intervals. For both, we take SEOBNRv2 
as the fiducial model because it has been calibrated against the 
highest number of high-accuracy NR simulations of aligned- 
spin binaries (38 in total, with mass ratio up to ~ 8), and is 
therefore likely to be the most accurate representation of true 
waveforms available at present). Overall, we find that (i) the 
two SEOBNR models (SEOBNRvl,2) disagree significantly 
for anti-aligned spinning binaries (matches below 80%), with 
their mismatches accumulating over lower frequency inspiral 
orbits; (ii) PhenomC and SEOBNRv2 produce drastically dif¬ 
ferent waveforms over most of the NSBH parameter space, 
except for the small mass-ratio + small spin corner; and (iii) 
both PN models show slightly better agreement with SEOB- 
NRv2 than PhenomC, but still restricted to small mass-ratios 
and small component spins, which is consistent with [37]. 


A. Faithfulness of models 

In Fig. 3, we examine the faithfulness of the two SEOBNR 
models. Neither of these models were used as templates in 
past LIGO searches, because they were published after ini¬ 
tial LIGO searches were completed, but both are promising 
candidate models for aLIGO. Focusing on stellar-mass NSBH 
binaries, we fix the NS mass to 1.4M©, the NS spin to 0, and 
allow the BH mass to vary over [3,15]M 0 and the BH spin 
to vary over the allowed range of SEOBNRvl [—1,0.6] [36]. 
We see that the agreement between the models is primarily 
influenced by the BH spin and secondarily by the mass ra- 
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Mass ratio Mass ratio 

FIG. 4. These plots are similar to Fig. 3, except they compare SEOBNRv2 with PhenomC (left panel) and TaylorF2 (right panel). From 
the left panel, we observe that PhenomC and SEOBNRv2 are faithful to each other for very small BH spins and q = mi/m 2 < 4. Their 
matches drop sharply for moderately spinning binaries and also above moderate mass-ratios. Also, the fall in matches of PhenomC with 
increasing mass-ratio is not monotonic. In the right panel, we find good agreement between TaylorF2 and SEOBNRv2 for binaries with very 
small BH spin. The best agreement, however, is for small anti-aligned spins, with the two models having < 95% matches for non-spinning 
binaries. Note that PhenomC uses TaylorF2 phasing prescription the early inspiral, including, however, only the leading and next-to-leading 
order spin-dependent terms. 



3456789 10 

Mass ratio 


FIG. 5. This figure is similar to Fig. 3, except it compares Tay- 
lorT4 with SEOBNRv2. We find that for mass ratios q > 3.5 or BH 
spin xbh ^ [—0 .1,0.6], the two models disagree significantly, with 
matches falling below 0.9, down to 0.4. 


tio. As expected, both agree in the comparable mass and 
non-spinning limits, where both incorporate information from 
NR simulations. We also find good agreement for aligned 
BH spins. However, when the BH spin is anti-aligned with 
the orbital angular momentum, SEOBNRvl and SEOBNRv2 
produce significantly different waveforms, with matches drop¬ 
ping below 0.8 for xbh < —0.5. This demonstrates that the 
more recent SEOBNRv2 model incorporates different spin- 
dependent phasing terms. However, to make statements about 
the accuracies of either, we must analyze both using high- 
accuracy NR waveforms, a comparison we turn to in the next 
section. 

In Fig. 4, we show the faithfulness of the PhenomC (left 
panel) and TaylorF2 (right panel) models against SEOBNRv2. 
We notice that PhenomC agrees with SEOBNRv2 only for 
very mildly spinning binaries with |xbh| <0.1 and compa¬ 


rable mass ratios q < 4. Over the remainder of the parameter 
space, the matches between PhenomC and SEOBNRv2 were 
found to be low. This is somewhat surprising, since PhenomC 
has been calibrated against spinning NR simulations with q 
up to 4 as well, albeit spanning fewer orbits and produced 
using a different NR code [72, 73]. Comparing with the right 
panel of Fig. 4, we find that TaylorF2 agrees with SEOBNRv2 
more closely for rapidly spinning NSBH systems with q < 3, 
as well as for binaries with small BH spin magnitude at all 
mass ratios. Since the inspiral portion of PhenomC phasing 
is the same as TaylorF2 (with the caveat that the former does 
not include the recently published spin-dependent 3 PN and 
3.5 PN contributions [65, 66]), we conclude that their differ¬ 
ences arise from the post-merger phasing prescription of Phe¬ 
nomC. We study this further later in this section, where we 
highlight the phase of binary coalescence where different ap- 
proximants disagree. 

We further show the faithfulness between TaylorT4 and 
SEOBNRv2 models in Fig. 5. We find that their faithfulness 
drops sharply with increasing mass ratio, falling below 0.9 for 
mass-ratios q > 4 for all values of BH spin. Only for near¬ 
equal mass low-spin binaries does TaylorT4 agree well with 
SEOBNRv2, which is consistent with past comparisons with 
NR simulations (e.g., Fig. 7 of [74]). Comparing with the 
right panel of Fig. 4, we see that (a) TaylorF2 has better over¬ 
all agreement with SEOBNRv2, and (b) TaylorF2 agrees bet¬ 
ter for small anti-aligned spins and TaylorT4 for small aligned 
spins. These differences are expected to decrease in future PN 
approximants, as higher order PN terms become available. 


B. Accumulation of mismatches with frequency 

For the adiabatic early-inspiral phase where the binary is 
well separated and inspirals relatively slowly, all GW models 
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FIG. 6. In this figure, we show the accumulation of mismatch between the two SEOBNR models that we consider in this paper by plotting 
the match (shown by colors) as a function of the low-frequency cutoff in the match calculation (horizontal axes) and black-hole spins (vertical 
axes). We choose 3 sets of NSBH systems corresponding to q — {5, 7,10} and ttibh = {7M©, 9.8M©, 14M©} (from left to right panels, 
respectively). We vary the black-hole spin (vertical axes) over the validity range for SEOBNRvl [36], xbh G [—1, 0.6]. We fix the mass of the 
neutron star to 1.4M© and the spin of the neutron star to zero. The solid, dotted, dash-dotted and dashed lines are contours of the frequencies 
starting from which the binary merges after 5,10, 20, 30 orbits, respectively. For anti-aligned BH spins, the matches are high when we integrate 
over the last few tens of orbits, but they fall significantly as we include more orbits (lower frequencies) in the computation. For moderate aligned 
BH spins, we find dephasing in the late-inspiral which is compensated for by lower frequency orbits. For comparison, the longest Numerical 
Relativity simulation to which either of the two models considered here have been calibrated to spans about 33 orbits [53, 67]. 
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FIG. 7. These figures are similar to Fig. 6, with the only difference that we compare here the PhenomC and SEOBNRv2 models. We observe 
that the two models agree over the last few (< 10) orbits, but their matches drop sharply over earlier late-inspiral orbits. The inclusion of very 
low-frequency orbits in match calculations leads to an increase in matches for |xbh| <0.2. The pattern shown in these figures suggests that 
dephasing is accrued rapidly close to the point where the model switches from TaylorF2 to its pre-merger phasing prescription. 


considered here are based on PN results. As the binary tight¬ 
ens, the PN approximation becomes less accurate. In order to 
capture the late-inspiral/plunge and merger phases, the IMR 
models either use purely phenomenological prescriptions or 
re-sum truncated PN results to add terms at all unknown 
higher orders in a controlled way. In both approaches, the 
model is calibrated against NR simulations of binary merg¬ 
ers, and models which are more extensively calibrated tend 
to be more robust. But do the different models agree in the 
late-inspiral regime, where the PN approximation is not valid 
and where we have no NR simulations available? To answer 
this question, we study here the mutual disagreement between 
models over different phases of binary coalescence. 

First, we examine the mismatch accumulation between 
the two EOB models for three representative sets of NSBH 
masses. In Fig. 6, the three panels correspond to mass-ratios 


q = mi/m 2 = {5, 7,10} (left to right), with the NS mass 
fixed at 1.4M©. The color shows the match between SEOB¬ 
NRvl and SEOBNRv2 as a function of the lower frequency 
cutoff on the match integral (shown on the x axes), for differ¬ 
ent values of BH spin (shown on the y axes). The four green 
curves on each panel are contours of the frequencies that mark 
“N orbits to merger”, with N = 30,20,10,5, respectively 
from left to right. We first note that the two models agree well 
for non-spinning binaries. When the BH spin is aligned to the 
orbital momentum, we observe some dephasing for binaries 
with xbh > 0.3, that accumulates over the last 30 or so orbits. 
When integrated over the entire waveform, this dephasing gets 
compensated for by the lower freqeuency orbits. For anti¬ 
aligned BH spins, however, the agreement between the two 
models is significantly low during the early inspiral. Over the 
last 15 — 20 orbits the two models have matches > 0.99, but 
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FIG. 8. These figures are similar to Fig. 6 with two differences: (a) here we compare the TaylorF2 and SEOBNRv2 models, and (b) the quantity 
shown is the match between the models as a function of the upper-frequency cutoff, with the low-frequency cutoff fixed at 15 Hz. Similar tc 
the trend observed in Ref. [37] for SEOBNRvl, TaylorF2 agrees with the SEOBNRv2 model at low frequencies, but their agreement dropj 
significantly during and after late-inspiral. The matches drop starting at lower frequencies with increasing BH spin magnitudes. These patterns 
are consistent with the fact that PN results decrease in accuracy with increasing orbital frequencies, especially as component spins becomes 
large. 
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FIG. 9. These figures are similar to Fig. 8 with the only difference that here we compare the TaylorT4 and SEOBNRv2 model. Comparing to 
Fig. 8, we find that (a) the patterns of mismatch accumulation in these figures are qualitatively similar to TaylorF2, and (b) the disagreement 
between TaylorT4 and SEOBNRv2 accumulates starting at lower frequencies and more drastically compared to TaylorF2. These results 
suggest that higher order PN terms in orbital phasing, especially spin dependent terms, are required to obtain a more accurate PN description 
of the late-inspiral phase. 


they drop below 0.95 around the 20 — 30 orbit mark, and fur¬ 
ther decrease monotonically as lower frequencies are included 
in the match calculation. Therefore, for anti-aligned spins, we 
would expect that NR simulations with lengths < 30 orbits, 
as used by Ref. [67] for SEOBNRv2, would match well with 
both the SEOBNRvl and SEOBNRv2 models, but that they 
would drastically disagree with at least one of the models 
early on. Given the lack of clear convergence of the SEOBNR 
models, we will investigate their accuracy in the following 
section by comparing them to long NR waveforms that probe 
the low-frequency regime where the models disagree. 

Next, we compare the PhenomC and SEOBNRv2 models 
by computing their matches with varying lower frequency cut¬ 
offs. The results are shown in Fig. 7, where all three panels 
are similar to Fig. 6 with the only difference that SEOBNRvl 
is replaced by PhenomC. We first note that the two models 
agree during the very early inspiral where PhenomC reduces 
to TaylorF2, as well as over the last few pre-merger orbits. 


Most of their dephasing accumulates in a relatively narrow 
frequency range centered at ^ 100 Hz, which is where Phe¬ 
nomC switches to its phenomenological phasing prescription. 
We also find that the spin on the BH affects the model agree¬ 
ment in two ways: (i) their dephasing increases with spin mag¬ 
nitude, and (ii) it also increases as the BH spin gets increas¬ 
ingly anti-aligned with the orbital angular momentum. 

Inspiral-only PN models have been shown to agree with 
SEOBNRvl during early inspiral and to monotonically di¬ 
verge as the orbital frequency rises [37]. To quantify their 
agreement with our fiducial model, SEOBNRv2, we compute 
matches between PN and SEOBNRv2 as a function of the up¬ 
per frequency cutoff (with the lower cutoff fixed at 15 Hz). 
In Fig. 8 we show the results for TaylorF2, with the upper 
frequency cutoff on the x-axes, BH spin on the y-axes and 
colors showing matches. The three panels correspond to the 
same representative NSBH systems as in Fig. 6. We find that 
for mildly anti-aligned BH spins with %bh > —0.2, Tay- 
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lorF2 agrees with SEOBNRv2 through most of late-inspiral. 
For other BH spin values, the two models start disagreeing 
at relatively low frequencies, e.g., for a q m 5 binary with 
Xbh = —0.6, the match drops to 0.9 between 15 Hz and 
300 Hz. In Fig. 9 we show the same results for TaylorT4. We 
find that for mild aligned BH spins with xbh < +0.2, Tay- 
lorT4 agrees with SEOBNRv2 through a significant portion of 
the inspiral. For higher or anti-aligned BH spins, the matches 
fall sharply below 0.7 as more of high frequency orbits are 
integrated over. The agreement gets restricted to even lower 
frequencies as the mass-ratio increases, for the entire range 
of BH spins. From this systematic frequency-dependent dis¬ 
crepancy between PN and SEOBNRv2, we expect that higher 
order spin-dependent PN corrections to orbital phasing are re¬ 
quired for a better PN modeling of the late-inspiral waveform. 

Finally, we note that some of the results presented in this 
section are qualitatively similar to [37], with the differences 
that (a) we additionally include the recently published spin- 
orbit tail (3PN) [65] and the next-to-next-to-leading order 
spin-orbit (3.5PN) contributions [66] in both of the PN mod¬ 
els, and (b) we include the SEOBNRv2 and PhenomC models 
here, both of which are capable of modeling binaries with very 
high black hole spins xbh — +1. 


VI. COMPARISON WITH NUMERICAL RELATIVITY 
SIMULATIONS 

In this section, we use our long NR simulations, described 
in Sec. II, to assess the accuracy of the different model wave¬ 
forms over the inspiral and merger phases of NSBH binary 
coalescences. 


A. Mismatch Accumulation for Models 

In order to quantify GW model errors over different phases 
of binary inspiral and merger, we compute matches between 
each model and NR waveforms over cumulative frequency in¬ 
tervals. The results are shown in Fig. 10- 14. Our main re¬ 
sults are as follows: (i) Of the two SEOBNR models, SEOB- 
NRv2 reproduces the late-inspiral and merger phases well for 
NSBH binaries with —0.9 < xbh < +0.6. In contrast, the 
SEOBNRvl model has an erroneous phase evolution during 
the late-inspiral phase, causing it to disagree both with the 
NR simulations and with SEOBNRv2 (cf. Fig. 6) (ii) The pre¬ 
merger phasing prescription of PhenomC does not reproduce 
NR waveforms well, as is confirmed by Ref. [75]. (iii) Of the 
two PN models we consider here, we found that TaylorT4 is 
more accurate for aligned BH spins, while TaylorF2 is more 
accurate for anti-aligned BH spins. 

In Figs. 10 and 11, we show the mismatches be¬ 
tween the TaylorT4, TaylorF2, PhenomC, SEOBNRvl, 
and SEOBNRv2 models and our aligned-spin simulations 
with q = rriBu/m^s = 9.8M 0 /1.4M 0 = 7, and 

XBH = {+0.6,+0.4} (ID SXS:BBH:202, SKS:BBH:203 
and SXS:BBH:204, c.f. table I), as a function of the lower 
(left panels) and upper (right panels) frequency cutoffs in the 


match calculation. First, we observe that SEOBNRv2 shows 
good agreement with NR with mismatches below 0.5% over 
all 55 — 88 orbits. We also find that SEOBNRvl agrees with 
NR over most of the inspiral orbits, but diverges closer to 
merger, with mismatches reaching 10% when considering the 
last few orbits. Therefore, we conclude that the disagreement 
between the two SEOBNR models for positive aligned spins 
that was seen in Fig. 6 stems from the phasing errors of SEOB¬ 
NRvl. Next, we find that the PhenomC model agrees well 
with the first few tens and the last few orbits of the NR wave¬ 
forms only. It accumulates significant phase errors in a narrow 
frequency band around 130 — 150 Hz, with mismatches rising 
above 10%. This explains the pattern of mismatch accumu¬ 
lation we observed in the middle panel of Fig. 7 and reaf¬ 
firms our conclusion that the pre-merger phasing prescription 
of PhenomC needs to be revisited for NSBH parameters [75]. 
We also find that both TaylorF2 and TaylorT4 (c.f. right pan¬ 
els of Fig. 10 and 11) agree with NR well up to the last 15 
or so pre-merger orbits. Closer to merger, their mismatches 
smoothly rise to 10%, which is expected as the PN approxima¬ 
tion degrades with increasing binary velocity. In addition, we 
find that TaylorT4 agrees with these NR waveforms to higher 
frequencies than TaylorF2, which is consistent with Fig. 5, 
where we show that TaylorT4 best agrees with SEOBNRv2 
for aligned, moderate BH spins. Finally, we note that the bot¬ 
tom row of Fig. 11 shows the comparison of different approxi- 
mants with a shorter NR simulation (ID SXS:BBH:203, start¬ 
ing frequency ^ 80 Hz) with the same physical parameters 
as a longer simulation (ID SXS:BBH:204, starting frequency 
~ 60 Hz). Therefore, the results for this simulation confirm 
those shown for the longer simulation (same figure, top row) 
for frequencies above 80 Hz. 

Next, we consider the case of NSBH binaries with anti¬ 
aligned BH spins. We show the mismatches between 
approximants and NR waveforms corresponding to q = 
^bhMns = 9.8M 0 /1.4M 0 = 7 andxBH = {— 1 0.4, —0.6} 
(ID SXS:BBH:205-207), in Fig. 12 and Fig. 13, respectively. 
We find that SEOBNRv2 shows the best agreement with NR, 
with mismatches below 0.2% across all orbits. SEOBNRvl 
monotonically accumulates increasing mismatches when we 
lower the lower frequency cutoff and appears to do the same 
over the later orbits when we increase the upper frequency cut¬ 
off, suggesting a phasing mismatch between the inspiral and 
merger portions of the waveform. Therefore, we conclude that 
the disagreement between the two SEOBNR models for anti¬ 
aligned spins that was seen in Fig. 6 stems from the phasing 
errors of SEOBNRvl. Similar to the aligned-spin cases, the 
PhenomC model agrees well with NR over the first few tens 
and the last couple of orbits, with mismatches below 1%, but 
accumulates large mismatches (10%) over a relatively narrow 
frequency range around 110 Hz. All of these observations are 
qualitatively similar to the aligned-spin cases. Lastly, we find 
that TaylorF2 shows excellent agreement with our NR wave¬ 
forms over the first 60 — 65 orbits, with mismatches below 
1%. It gradually diverges during the last 5 pre-merger orbits, 
with high mismatches, which is consistent with the right panel 
of Fig. 4. TaylorT4, on the other hand, performs worse, with 
mismatches accumulating earlier on and rising to 10%. The 
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FIG. 10. For a NSBH binary with q = mBu/m^s = 9.8M©/1.4M© = 7 and xbh = +0.6, these figures show the mismatch of TaylorF2, 
TaylorT4, SEOBNRvl, SEOBNRv2 and PhenomC waveforms against our simulation ID SXS:BBH:202 (see Table I) as a function of the lower 
(left) and upper (right) frequency cutoff on the overlap integral. The NR waveform starts from GW frequency /o — 80 Hz. The dashed curve 
with dark grey shading (in both panels) shows mismatches because of the tapering and high-pass filtering of NR waveforms, which we do to 
reduce Gibbs phenomena and spectral leakage upon Fourier transformation. Because of the width of the tapering window, we begin filtering at 
82.5 Hz. The solid curve with light grey shading (only barely visible at the bottom of the left panel and below the range of mismatches shown 
in the right panel) shows the mismatches between NR waveforms at the highest and second-highest available numerical resolutions. 
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FIG. 11. The top two panels of this figure are similar to the two panels of Fig. 10, with the difference that the system considered here 
has q = ttibh/jtjns = 9.8M©/1.4M© = 7, xbh = +0.4, and starts at GW frequency ~ 60 Hz. This corresponds to simulation 
ID SXS:BBH:204 (see Table I). The bottom two panels are similar to the top two with the difference that these correspond to simulation 
ID SXS:BBH:203, which has the same physical parameters as ID SXS:BBH:203 but a higher starting GW frequency, i.e. ~ 80 Hz. Because 
of the width of the tapering windows, we begin filtering at 63.5 Hz and 83 Hz, respectively. 







































































































12 


Orbits to merger 

71.9 49.6 32.4 20.4 12.9 7.8 4.5 2.6 1.5 1.0 



Orbits to merger 

£1.9 49.6 32.4 20.4 12.9 7.8 4.5 2.6 1.5 1.0 



60 72 88 108 131 159 194 236 287 350 


Lower Frequency cutoff (Hz) 
Orbits to merger 



Upper Frequency cutoff (Hz) 

Orbits to merger 

£1.9 49.6 32.4 20.4 12.9 7.8 4.5 2.6 1.5 1.0 



60 72 88 108 131 159 194 236 287 350 


Lower Frequency cutoff (Hz) 


Upper Frequency cutoff (Hz) 


FIG. 12. The top two panels of this figure are similar to the two panels of Fig. 10, with the difference that the system considered here 
has q m mBu/m^s = 9.8M©/1.4M© as 7, xbh = —0.4, and starts at GW frequency ~ 60 Hz. This corresponds to simulation 
ID SXS:BBH:206 (see Table I). The bottom two panels are similar to the top two with the difference that these correspond to simulation 
ID SXS:BBH:205, which has the same physical parameters as ID SXS:BBH:206 but a higher starting GW frequency, i.e. ~ 80 Hz. Because 
of the width of the tapering windows, we begin filtering at 63.5 Hz and 84 Hz, respectively. 


Orbits to merger 



Orbits to merger 
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FIG. 13. These figures are similar to the two panels of Fig. 10, with the difference that the system considered here has q = ttz-bh/ht-ns = 
9.8M©/1.4M© = 7, xbh = —0.6, and starts at GW frequency ~ 80 Hz. This corresponds to simulation ID SXS:BBH:207 (see Table I). 
Because of the width of the tapering window, we begin filtering at 88 Hz. 
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Orbits to merger 



Orbits to merger 
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FIG. 14. These figures are similar to the two panels of Fig. 10, with the difference that the system considered here has q = mBa/m^s — 
7M©/1.4M© = 5, xbh = —0.9, and starts at GW frequency ~ 80 Hz. This corresponds to simulation ID SXS:BBH:208 (see Table I). 
Because of the width of the tapering window, we begin filtering at 86 Hz. 
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FIG. 15. The top panels in this figure are reproductions of the top two panels of Fig. 11, and are shown here for direct comparison. The bottom 
two panels show identically computed quantities, with the only difference from the top two being that the NR waveform used corresponds 
to the second-highest numerical resolution instead of the highest. The case shown here has the highest mismatches between the highest and 
second-highest resolution NR waveforms, and therefore serves as a conservative example of the robustness of our results to NR errors. 


bottom row of Fig. 12 shows identical comparisons with the 
shorter NR waveform ID SXS:BBH:205 with the same phys¬ 
ical parameters as ID SXS:BBH:206. As in the aligned case, 
these panels agree with and provide a confirmation for the re¬ 
sults shown in the top row of the figure, for frequencies above 
80 Hz. 

Finally, we compare different approximants with NR sim¬ 
ulation ID SXS:BBH:208, which has a smaller mass-ratio 
q — 5 but the most extremal BH spin of all the cases we 
consider, with xbh = —0.9. Compared to the above two anti¬ 
aligned spin cases, we find that the SEOBNRvl mismatches 
rise to 5%, indicating that the phasing errors of SEOBNRvl 


get worse with BH spin magnitude. SEOBNRv2 still repro¬ 
duces the NR waveform the best, while PhenomC and both PN 
approximants show similar patterns to the q = 7 anti-aligned 
spin cases. 

Overall, we conclude that the more recently calibrated 
SEOBNRv2 [67] model reproduces the late-inspiral and 
merger phases well for NSBH binaries with —0.9 < xbh < 
+0.6. This model was calibrated to 8 non-spinning and 30 
non-precessing spinning simulations with mass-ratios up to 
q = 8, and most of the NSBH systems we consider here are 
within the calibration range of the model. Therefore, we con¬ 
clude that SEOBNR performs well when interpolated within 
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its calibration range. We also conclude that PhenomC does not 
reproduce NR waveforms well, and its NR-calibrated portion 
needs to be revisited [75]. Of the two PN models we con¬ 
sider here, TaylorT4 and TaylorF2, we found that the former 
has better accuracy for aligned-spin binaries, while the latter 
is more accurate for the case of anti-aligned BH spins. These 
patterns, however, are likely coincidental. 

Finally, note that in detection searches there is an additional 
degree of freedom corresponding to the maximization of the 
SNR over the intrinsic waveform parameters. We investigate 
the suitability of the GW models considered here as detection 
templates in Sec. VII. 


B. Impact of numerical errors in NR simulations 

Fig. 2 shows that some of the NR simulations (especially 
the ones with %bh = ±0.4) do not show explicit conver¬ 
gence with an increase in numerical grid resolution. There¬ 
fore, one might question the effect of the errors in the nu¬ 
merical waveforms that arise because of a finite grid reso¬ 
lution. In order to quantify this, in this section, we repeat 
the overlap calculations for the xbh = +0.4 configuration, 
using the NR waveform produced at the second-highest grid 
resolution. We choose this configuration because it exhibits 
the highest NR(highest resolution)-NR(next-to-highest reso¬ 
lution) mismatches, as shown by the solid line bounding the 
light grey region in Fig. 11 . We show the results in the bottom 
row of Fig. 15. The top row in the figure is a reproduction 
of the top row of Fig. 11, which uses the NR waveform pro¬ 
duced at the highest grid resolution. From the left two pan¬ 
els, we first notice that the PhenomC, TaylorT4 and TaylorF2 
mismatches are sufficiently large for the variations because 
of NR waveform differences to be inconsequent. Further, for 
both SEOBNRvl and SEOBNRv2, we find a small change 
in the mismatches near the left and the right edge of the two 
left panels. These fluctuations are entirely consistent with the 
value of the NR-NR mismatch (shown by the solid grey line), 
when we apply the triangle inequality to the square roots of 
the SEOBNR-NR and NR-NR mismatches (see Sec. Ill of 
Ref. [76] for a brief discussion on manipulation of waveform 
mismatches arising from independent sources). The same is 
true when we compare the two right panels in Fig. 15. In this 
case, the changes visible in the SEOBNR-NR mismatches are 
below the upper bound set by adding the square roots of the 
SEOBNR-NR and NR-NR mismatches. 

We therefore conclude that the fluctuations in the mis¬ 
matches computed in the previous subsection are within the 
error bounds set by comparing NR waveforms at the highest 
and second-highest numerical grid resolutions. These bounds 
are shown explicitly in light grey shading bounded by solid 
lines in all panels of Fig. 10- 13. As these fluctuations remain 
below 0.3% in all cases, our conclusions remain robust to NR 
waveform errors. 



FIG. 16. In this figure, we show the ineffectualness of different wave¬ 
form approximants against our NR waveforms, as a function of black 
hole spin. Dashed lines join points corresponding to the cases where 
the NR waveform starts at ~ 80 Hz (ID: 01, 02a, 03a, 04), while the 
solid lines correspond to the cases where the NR waveform starts at 
~ 60 Hz (ID: 02b, 03b). The isolated points correspond to the q = 5 
simulation (ID: 5). We find that SEOBNRv2 consistently recovers 
the highest fraction (> 99.8%) of the optimal SNR when optimizing 
over intrinsic binary parameters. SEOBNRvl is also fairly effec¬ 
tual. The Taylor approximants recover 97 — 98% of the SNR and 
noticeably more for the longer NR waveforms (for fixed binary pa¬ 
rameters). This is the expected trend, as longer waveforms extend to 
lower frequencies, where the PN approximation is better. PhenomC 
is effectual against NR for xbh > 0, recovering > 99% of the SNR. 
For anti-aligned spins, however, its effectualness was found to be low 
and decreasing with increasing NR waveform length. 

VII. EFFECTUALNESS AND PARAMETER BIAS 


While accurate parameter estimation requires that wave¬ 
form models be faithful to the true signal for given binary pa¬ 
rameters, detection searches allow for the additional degree of 
freedom of maximizing the SNR over intrinsic binary param¬ 
eters. With this freedom, intrinsic waveform model errors can 
be compensated by slight shifts in the physical parameters. As 
a result, a GW signal will be better matched with a template 
waveform with slightly incorrect physical parameters. In this 
section, we investigate the recovery of SNR using different 
approximants and the associated biases in the recovered max¬ 
imum likelihood estimates for binary masses and spins. We 
take all of our NR waveforms, produced at the highest numer¬ 
ical resolution, and compute their overlaps against 500,000 
modeled waveforms with physical parameters in the vicinity 
of the true NR parameters. We compute waveform overlaps 
integrating from the starting frequency of the simulation in 
question, up to the Nyquist frequency corresponding to the 
waveform sampling rate. 
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FIG. 17. These figures show the fractional error in the recovered maximum likelihood parameters as a function of black hole spin, except for 
black hole spin for which the actual error value is shown. Each approximant is denoted with a unique color. As in Fig. 16, dashed lines join 
points corresponding to NR waveforms that start at ~ 80 Hz (ID: 01, 02a, 03a, 04), while the solid lines correspond to the ones starting at 
~ 60 Hz (ID: 02b, 03b). The isolated points correspond to the q = 5 simulation (ID: 5). While the left-to-right trends show us the effect of 
component spin, comparing dashed and solid lines show us the effect of including more inspiral cycles. 


In Fig. 16 we show the fractional loss in SNR, or the inef¬ 
fectualness of approximants, as a function of BH spin. The 
q = 7 cases are connected with straight lines, while the single 
q = 5 case is shown by a point. Each approximant is shown 
with a different color. Solid and dashed lines join points cor¬ 
responding to ID 2a, 3a and ID 1,26,36,4, respectively. De¬ 
tection searches use banks of templates that correspond to a 
grid in the parameter space, and the SNR loss due to the dis¬ 
creteness of the grid will be in addition to those shown in 
this figure. We find that both EOB models recover more than 
99.5% of the optimal SNR for both aligned and anti-aligned 
BH spins, despite different faithfulness. This is a good exam¬ 
ple of a shift in waveform parameters compensating for model 
phasing errors. For aligned-spin cases, both PN approximants 
as well as PhenomC are effectual with SNR recovery above 
98%. For all of the anti-aligned spin cases, however, we found 
relatively low SNR recovery using PN approximants, which is 
to be expected as anti-aligned spin binaries merge at lower fre¬ 
quencies than aligned-spin cases and therefore have relatively 
less signal power in the inspiral cycles above a given physi¬ 
cal frequency (here 60 — 80 Hz). PhenomC also shows sig¬ 
nificantly low effectualness for anti-aligned systems, recover¬ 
ing < 96% of the SNR for q = 7 cases and 91.5% for the 
q = 5, xbh = —0.9 case. We therefore conclude that (a) both 
SEOBNR models are sufficiently accurate to model NSBH 
templates in aLIGO detection searches, and (b) PN approxi¬ 
mants are also viable for use as filter templates for aligned- 


spin NSBH systems. 

In Fig. 17 we show the fractional difference between the 
parameters that maximize the SNR recovery for each of the 
approximants, i.e. the maximum likelihood parameters, and 
the true physical parameters of the system. Table II lists the 
same for the IMR approximants, and Table III for the inspiral- 
only approximants. These differences quantify the systematic 
bias intrinsic to each approximant. From the top left panel of 
Fig. 17, we observe that as the number of waveform cycles in¬ 
creases (monotonically with BH spin), so does the accuracy of 
the recovered chirp mass. For aligned spins, all approximants 
but SEOBNRv2 converge at a systematic —1% bias in chirp 
mass recovery. SEOBNRv2 shows a smaller (< 0.8%) bias in 
chirp mass for all BH spins. In the top center panel of Fig. 17, 
we show the bias in mass-ratio recovery. The spin-mass-ratio 
degeneracy that enters at the sub-leading order in PN phasing 
is manifest here, and for PN models, we find i) a larger bias in 
rj for aligned spins that increases with the spin magnitude and 
ii) a smaller (or negative) bias for anti-aligned spins. We see 
a similar trend for the IMR approximants, with SEOBNRv2 
showing the smallest systematic bias. We also show biases 
in other binary mass combinations, total mass in the top right 
panel, BH mass in the bottom left, and NS mass in the bot¬ 
tom center. These show that none of the individual masses or 
their sum are nearly as accurately measured as the chirp mass 
of the binary. Additionally, the mass of the smaller compo¬ 
nent, here the NS, is slightly more biased than for the more 
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SEOBNRvl 


Actual Sbh 

M c 

r] m i 

m 2 

Sbh 

0.6 

-1.12% 

50.45% -29.9% 

28.6% 

-0.184 

0.4 (80Hz) 

-1.02% 

29.1% -19.5% 

15.6% 

-0.180 

0.4 (60Hz) 

-0.82% 

72.1% -38.8% 

44.1% 

-0.355 

-0.4 (80Hz) 

0.16% 

6.19% -4.40% 

3.67% 

-0.066 

-0.4 (60Hz) 

-0.44% 

31.3% -20.2% 

17.6% 

-0.382 

-0.6 

0.97% 

-16.1% 15.1% 

-8.29% 

0.234 

-0.9 (q = 5) 

0 .66% 

-5.2% 5.3% 

-2.7% 

0.102 

SEOBNRv2 

Actual Sbh 

Me 

r] m\ 

7712 

Sbh 

0.6 

-0.49% 

36.3% -22.8% 

20.5% 

-0.0966 

0.4 (80Hz) 

-0.16% 

-3.48% 2.57% 

-2.13% 

0.00808 

0.4 (60Hz) 

-0.10% 

0.54% -0.51% 

0 .21% 

0.00444 

-0.4 (80Hz) 

0 .10% 

-7.84% 6.48% 

-4.35% 

0.0836 

-0.4 (60Hz) 

-0.18% 

6.64% -5.03% 

3.58% 

-0.0864 

-0.6 

-0.77% 

13.5% -10.2% 

6.87% 

0.234 

-0.9 (q m 5) 

0.7% 

-22.7% 23.7% 

-13.7% 

0.387 

PhenomC 

Actual Sbh 

Me 

r] mi 

7772 

Sbb 

0.6 

-0.82% 

70.5% -38.1% 

42.8% 

-0.191 

0.4 (80Hz) 

-1.18% 

67.3% -37.0% 

40.0% 

-0.353 

0.4 (60Hz) 

-1.14% 

69.6% -38.0% 

41.7% 

-0.373 

-0.4 (80Hz) 

6.59% 

39.7% -19.1% 

31.3% 

-0.0259 

-0.4 (60Hz) 

5.61% 

19.1% -8.06% 

17.1% 

0.202 

-0.6 

9.90% 

45.9% -19.8% 

39.7% 

0.0366 

-0.9 (q = 5) 

9.6% 

-5.4% 14.8% 

5.8% 

0.693 


TABLE II. In this table we list the fractional error in the recovered 
maximum likelihood parameters for different IMR approximants. 
For BH spin, we give the actual difference between maximum likeli¬ 
hood and true parameter. 



FIG. 18. In this figure, we show the lowest SNR value below which 
a modeled waveform (using the respective approximant) with the 
same parameters as the true signal waveform will be indistinguish¬ 
able from the true signal waveform. Here, the true signal waveform 
is represented by the corresponding NR waveform. We use the cri¬ 
terion proposed in Ref. [49] for this calculation. As in Fig. 16, solid 
lines correspond to the two longer simulations that start at 60 Hz, 
dashed lines are for the shorter simulations that start at 80 Hz, and 
the isolated points correspond to the q = 5 simulation (ID: 5). For 
all but SEOBNRv2 (and marginally SEOBNRvl), this threshold is 
below the SNR cutoff (= 5.5, shown by horizontal black line) em¬ 
ployed by past LIGO-Virgo searches, and therefore their use would 
likely degrade the extraction of information from detector data. 


TaylorT4 


Actual Sbh 

M c 77 

mi 

7772 

Sbh 

0.6 

-0.67% 81.9% 

-42.4% 

52.1% 

-0.207 

0.4 (80Hz) 

-1.73% 50.0% 

-30.1% 

27.5% 

-0.366 

0.4 (60Hz) 

-1.23% 16.92% 

-12.7% 

8.31% 

-0.191 

-0.4 (80Hz) 

-4.35% -34.2% 

29.5% 

-23.2% 

-0.109 

-0.4 (60Hz) 

-3.95% -23.0% 

16.5% 

-16.6% 

-0.254 

- 0.6 

-4.09% -41.1% 

40.2% 

-27.1% 

-0.00606 

-0.9 (q m 5) 

-3.94% -36.5% 

36.6% 

-26% 

-0.0374 

TaylorF2 

Actual Sbh 

M C 77 

mi 

7772 

Sbh 

0.6 

-1.30% 79.1 

-41.8% 

48.9% 

-0.281 

0.4 (80Hz) 

-0.67% 79.6% 

-41.5% 

50.1% 

-0.398 

0.4 (60Hz) 

-0.82% 72.1% 

-38.8% 

44.1% 

-0.355 

-0.4 (80Hz) 

2.02% 5.61% 

-2.20% 

5.27% 

0.0828 

-0.4 (60Hz) 

1 .12% -8.82% 

8.42% 

-3.94% 

0.172 

- 0.6 

3.04% -4.53% 

6.74% 

0.40% 

0.247 

-0.9 (q = 5) 

2.16% -41.5% 

54% 

-24.6% 

0.735 


TABLE III. This table is identical to Table II, but for inspiral-only 
approximants. 


massive component (here the BH). Finally, the bottom right 
panel of the same figure shows the bias in recovered values of 
black hole spin. We observe a systematic underestimation of 
BH spin when it is aligned with the orbit, and smaller under¬ 
estimation or overestimation when it is anti-aligned. Both of 


these patterns are exacerbated with BH spin magnitude. 

Overall, we found SEOBNRv2 to have the smallest system¬ 
atic biases in parameter recovery. Therefore, we recommend 
its use in aLIGO parameter estimation efforts focused on non- 
precessing NSBH binaries. 

Having ascertained the systematic parameter biases that are 
applicable in the limit of high SNR, we ask the question: how 
loud does an incoming GW signal have to be before a mod¬ 
eled waveform with the same parameters is distinguishable 
from the true, measured waveform? Ref. [49] proposed the 
criterion: (Sh\Sh) < 1 , where Sh = h true — h a pprox , which 
is sufficient for proving the indistinguishability of the mod¬ 
eled waveform h a PP rox from the true signal h true . We use it 
to calculate the effective SNR p e ff below which different ap¬ 
proximants are indistinguishable from true (NR) waveforms, 
and show it in Fig. 18 as a function of black hole spin, for 
all <7 = 7 and q = 5 cases. We immediately observe that 
for TaylorF2, TaylorT4, and PhenomC, the SNR threshold for 
distinguishability is below what is chosen as the single detec¬ 
tor SNR lower cutoff in LIGO searches (= 5.5) and there¬ 
fore using these approximants would likely degrade scientific 
measurements. For SEOBNRvl, inclusion of more inspiral 
cycles for anti-aligned spin cases ID 3,5 cause a drop in p e # 
from 12 — 6. This is consistent with Fig. 12 which shows that 
SEOBNRvl monotonically diverges from the reference NR 
waveform(s) when more of inspiral cycles are considered. For 













































17 


aligned-spin cases, we find that SEOBNRvl is always distin¬ 
guishable from a real signal with SNR above the lower cutoff 
for aLIGO searches. SEOBNRv2 is indistinguishable from 
true waveforms to fairly high SNRs ~ 15 — 20 for anti-aligned 
spins, but this threshold lowers when we consider the longer, 
aligned-spin cases. 

Lastly, we note that the SNR in consideration here is in¬ 
tegrated from the starting frequency of the NR waveform in 
question, i.e. ^ 60 or 80 Hz, and corresponds to about 
50 — 70% of the total SNR accumulated over the entire aLIGO 
frequency band starting from 15 Hz. These results are there¬ 
fore likely to be pessimistic for PN models, since these models 
do better at lower frequencies, where PN is more accurate. 


VIII. CONCLUSIONS 

With the first observations of the Advanced LIGO and 
Virgo detectors imminent, rapid development of data anal¬ 
ysis methods is underway within the LIGO Scientific Col¬ 
laboration and the Virgo Collaboration. Gravitational-wave 
astronomers are shaping matched-filtering-based targeted 
searches for neutron star - black hole binaries. As a step 
forward from most of earlier LIGO-Virgo searches (which 
used non-spinning waveforms as templates, e.g. [24-27]), Ad¬ 
vanced LIGO searches plan to employ aligned-spin wave¬ 
forms as templates. This is motivated towards increasing the 
sensitivity of the searches for binaries with spinning compo¬ 
nents. Recent work has shown that even if the component 
spins are not aligned to the orbital angular momentum and 
result in orbital precession, aligned-spin waveform templates 
would likely have significantly better sensitivity towards such 
systems than non-spinning waveform templates [33]. 

Recent progress in numerical relativity has allowed for 
faster and more accurate general-relativistic numerical sim¬ 
ulations of inspiraling black holes, including the effect of 
component spin [54]. With these advances, more and longer 
simulations of compact binary motion have become possi¬ 
ble [77]. While the possibility of using numerical relativ¬ 
ity waveforms directly as search templates has been demon¬ 
strated [76], Bayesian parameter estimation efforts will re¬ 
quire the ability to generate template waveforms for arbi¬ 
trary source parameters. This is computationally prohibitive 
with the current NR technology, and therefore approximate 
waveform models are indispensible. Using strong-field infor¬ 
mation from numerical relativity, Effective-One-Body (EOB) 
and phenomenological (Phenom) models have been devel¬ 
oped and calibrated to accurately model the late-inspiral mo¬ 
tion of compact binaries all the way through merger. The NR 
input here has been critical, since the post-Newtonian expres¬ 
sions that form the basis of all IMR models are perturbative 
expansions in the invariant velocity v/c, which become inac¬ 
curate in the strong field, rapid motion regime. While the Phe¬ 
nom models are closed-form in frequency domain and there¬ 
fore the least expensive to generate, reduced-order methods 
have been recently applied to the EOB family to mitigate their 
computational cost [78]. 

In this paper, we present 7 new NR simulations, 6 with 


q = m 1 /m 2 = 7 and black hole spins xbh = {±0.4, ±0.6}, 
and 1 with q = 5 and xbh = —0.9. The spin of the smaller 
object (a black hole representing the neutron star) is held at 
0. For xbh = ±0.4, we perform two simulations each, 
one starting at a gravitational-wave frequency of 60 Hz and 
the other starting at 80 Hz (corresponding to a total mass of 
1.4M 0 ± 7 x 1.4M 0 = 11.2M 0 ). These span 36 — 88 pre¬ 
merger orbits. For all other parameter values, our simulations 
start close to 80 Hz when scaled to appropriate NSBH masses 
(cf. Table I). Using these simulations, we study the accuracy 
of different waveform approximants and their effectualness as 
models for search and parameter estimation templates for the 
Advanced detector era. 

Our investigation of the faithfulness of two inspiral-merger- 
ringdown models and two inspiral-only PN models shows that 
both PN models become increasingly unfaithful with increas¬ 
ing BH spin magnitudes as well as with binary mass ratio, 
with overlaps falling below 50%. This is consistent with 
a similar study [37] and is indicative of the breakdown of 
the PN approximation with increasing binary velocity. We 
find that PhenomC disagrees in an even larger portion of the 
parameter space with SEOBNRv2, with overlaps above 0.9 
for near-equal mass binaries with spin magnitude below 0.3. 
Somewhat surprisingly, we also find that the two SEOBNR 
models diverge significantly for anti-aligned BH spins. Next, 
we investigate the GW frequency dependence of these model 
disagreements to disambiguate the portion of the binary coa¬ 
lescence process, that each model fails to capture accurately. 
As expected, the PN models describe the inspiral well but 
break down closer to merger. We found that PhenomC ac¬ 
cumulates most of the mismatch against SEOBNRv2 close 
to the frequencies where it switches its phase prescription 
from one piece of a piece-wise continuous function to another. 
Lastly, for aligned spins, SEOBNRvl accumulates phase dif¬ 
ferences close to merger, but for anti-aligned spins it agrees 
with SEOBNRv2 close to merger, with most of its mismatch 
being accumulated earlier on during the late inspiral. We 
present these results in detail in Sec. V. 

When we study the mismatch accumulation of GW mod¬ 
els as a function of frequency, we find that the PN models are 
faithful at the lower frequencies of our NR waveforms, and di¬ 
verge close to merger. We find that PhenomC reproduces NR 
waveforms accurately during the last 2 — 5 and first 20 — 30 
orbits, but accumulates significant mismatches over a span of 
~ 20 orbits around 100 Hz. For SEOBNRvl, we find that 
(a) for aligned-spin binaries, it slowly accumulates phase er¬ 
ror over the last ~ 5 orbits with mismatches rising to 10%, 
but agrees well earlier on, but (b) for anti-aligned binaries, it 
agrees well with NR during the last 10 — 20 orbits and di¬ 
verges monotonically when more inspiral orbits are included, 
with mismatches rising to 2%. In contrast, SEOBNRv2 repro¬ 
duced all of our NR waveforms well, throughout the probed 
frequency range, with mismatches below 1%. We summarize 
these results in Sec. VI. 

Further, we study the effectualness of these models as de¬ 
tection templates. Detection searches allow for the additional 
freedom of maximizing the match of template with signal over 
the intrinsic parameters of the templates, allowing for partial 
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compensation of modeling errors. Using this freedom, we 
computed the mass-and-spin optimized overlaps between our 
set of numerical waveforms and those generated using differ¬ 
ent waveform approximants. For BH spins aligned with the 
orbital momentum, all models are effectual against our NR 
waveforms, with SNR recovery above 98%. For anti-aligned 
BH spins, we find that both PN models are less effectual, 
with SNR recovery dropping below < 96.5%. This is ex¬ 
pected, since our NR waveforms have fewer inspiral orbits for 
anti-aligned spins, even though they span similar frequency 
ranges. We find that PhenomC recovers the lowest SNRs (be¬ 
low 90 — 95%, depending on spin) for anti-aligned spins. We 
find that both SEOBNR models recover more than 99.8% of 
the optimal SNR when maximized over physical parameters. 
Therefore, we recommend using SEOBNR models to model 
non-precessing templates in aLIGO detection searches. We 
also note that the unpublished PhenomD model (an improve¬ 
ment over PhenomC) has shown promising results in terms 
of accurately capturing the merger waveforms for high mass- 
ratio non-precessing binaries [79], and would therefore be an¬ 
other suitable candidate for modeling search and parameter 
estimation templates in aLIGO era. 

Finally, we also investigate the systematic bias in maxi¬ 
mum likelihood recovered parameters that GW models will 
incur if used to model parameter estimation templates. We 
find that while the chirp mass is recovered increasingly and 
very accurately (within a percent) with increasing number of 
binary orbits in detector frequency band, the spin-mass-ratio 
degeneracy makes accurate determination of other parameters 
more difficult. For mass-ratio, we find that TaylorT4, Tay- 
lorF2, PhenomC, and SEOBNRvl have a 20 — 50 + % bias, 
with SEOBNRv2 relatively the most faithful to NR with a 
2 — 38% bias. We also find that the models consistently under¬ 
estimate BH spins for aligned-spin binaries, and overestimate 
(or slightly underestimate) it for anti-aligned binaries. Over¬ 
all, we find SEOBNRv2 to be the most faithful approximant 
to NR. As these biases are applicable in the high SNR limit, 
we also investigate the SNR limit below which different ap¬ 
proximants are essentially indistinguishable from NR wave¬ 
forms [49]. We find that for PN and Phenom models, this is 
never the case for signals with SNRs above 5; but up to SNR 
^ 10 — 22 (depending on spin), any further increase in 
the accuracy of SEOBNRv2 will not affect the extraction of 
scientific information from detector data. We describe these 
results in Sec. VII. Therefore, for NSBH binaries with mod¬ 
erate spins, parameter estimation efforts will benefit from us¬ 
ing SEOBNRv2 templates in aLIGO era. However, given the 
drop in accuracy of the SEOBNRvl model outside its range of 
calibration, we also recommend further investigating SEOB- 
NRv2 at more extremal component spins, and we recommend 
trusting SEOBNRv2 only within its calibration range. 

We note that we ignore a good fraction (35 — 45%, depend¬ 
ing on BH spin) of the signal power by considering only fre¬ 
quencies above 60 — 80 Hz. Therefore, our results are accu¬ 
rate in the high-frequency limit and are likely to be pessimistic 


for the PN-based inspiral-only models. We plan to extended 
the study to lower frequencies in future work. We also note 
that there is tremendous ongoing effort to model the effect 
of the tidal distortion of neutron stars during NS-BH merg¬ 
ers [46]. They are expected to be measurable with aLIGO 
detectors [45]. However, matter effects affect binary phasing 
at 5+PN order, while there are lower order unknown spin de¬ 
pendent terms in PN phasing whose lack of knowledge will 
have a larger impact on the detection problem [80] . This mo¬ 
tivates our choice to ignore neutron star matter effects in our 
numerical simulations, and instead treating the neutron stars 
as low-mass black holes. We have also ignored the effect of 
the possible tidal disruption of NS in this study, and leave its 
detailed analysis to future work. However, since NS disrup¬ 
tion affects NSBH waveforms only at very high (>1.2 kHz) 
frequencies [47] where aLIGO has significantly reduced sen¬ 
sitivity, ignoring it is unlikely to affect the accuracy of our 
analysis. Finally, we note that we consider only the dominant 
l = =bm = 2 waveform multipoles in this study, since (i) 
other multipoles have much smaller (by orders of magnitude) 
contribution to the SNR, and (ii) none of the IMR models con¬ 
sidered here include sub-dominant modes. 
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